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A compressible CFD code has been developed to study the heat transfer characteristics of a /3-type Stirling 
engine with a very simple design and geometry. The results include temperature contours, velocity vec¬ 
tors, and distributions of local heat flux along solid boundaries at several important time steps as well as 
variations of average temperatures, integrated rates of heat input, heat output, and engine power. It is 
found that impingement is the major heat transfer mechanism in the expansion and compression 
chamber, and the temperature distribution is highly non-uniform across the engine at any given moment. 
The results, especially the rates of heat transfer, are quite different from those obtained by a second-order 
model. The variations of heat transfer rates are much more complicated than the simple variations 
returned by the second-order model. This study sheds light into the complex heat transfer mechanism 
inside the Stirling engine and is very helpful to the understanding of the fundamental process of the 
engine cycle. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

Stirling engine, a nearly 200 years old invention, has attracted 
much attention and been the subject of intensive research by many 
industrial and academic institutes alike in recent decades due to its 
many specific advantages [1-4], These advantages are: very high 
efficiency, silence in operation, low level of toxic emission if 
powered by fuel, and the ability to use almost any kind of heat 
sources, including solar, biological, geothermal, or even industrial 
waste heat. Nowadays, it has been regarded as an important green 
energy device and potentially one of the most promising solutions 
to the global warming problem. Among its most successful applica¬ 
tions is the micro combined heat and power (CHP) systems where 
the power output of the engine is used to generate electricity, and 
the waste heat from engine cycle is used for heating purpose [5-7], 
However, there are also disadvantages associated with Stirling 
engines. It is low in the power/weight ratio, unable to change tor¬ 
que and power swiftly, and expensive (especially for those using 
concentration solar energy); and these have precluded the engine’s 
spread to more general applications. Despite these disadvantages, 
Stirling engine technology has been constantly improved to 
produce new generation of engines with higher power/weight 
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ratio, efficiency, and power output; and Stirling engines might 
become more popular and spread to more general applications in 
the future. 

Improvement on Stirling engine’s performance relies heavily on 
extensive numerical and experimental research. To design a new 
Stirling engine, a numerical analysis often precedes experiment 
to find some general directions to follow. It is also a more suitable 
tool than experiment to carry out parameter optimization. There¬ 
fore, many numerical models for Stirling engine cycle have been 
developed. A detailed overview of the existing numerical models 
for Stirling engine cycles can be found in Mahkamov [8]. These 
models are broadly classified as the first-, second-, and third-order 
models, according to the model hierarchy. Here, the term “order” 
has nothing to do with the mathematical terms in the governing 
equations. In these models, the entire engine is divided into 3, or 
5, or more sections (control volumes), so they are practically 
zero- or one-dimensional models. Among them, the model pro¬ 
posed by Schmidt is the stereotype, and improved variants of 
Schmidt’s model are nowadays the most widely adopted Stirling 
engine models [9-11 j. However, because these models are zero- 
or one-dimensional models, they, at best, only resolve spatial vari¬ 
ations in the axial direction of the engine. In this case, any variation 
in the transverse direction that is important for predicting the 
engine’s performance accurately will not be accounted for. Neither 
can they resolve the effects caused by changes in some geometrical 
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Nomenclature 



Cp 

constant pressure specific heat (kj kg -1 K -1 ) 

t 

time (s) 

G 

regenerative channel gap (m) 

tp 

time period of an engine cycle 

h 

heat transfer coefficient (W K -1 ) 

T 

temperature (K) 

k 

thermal conductivity (W nr 1 I< -1 ) 

U, V 

velocity components respectively in x- and r-direc- 

la 

height of displacer (m) 


tions (ms -1 ) 

U 

height of the engine (m) 

u^v 

local frame moving velocity components respectively 

hi h. 

1 3 , U lengths of the linkages of the rhombic drive mecha- 


in x- and r-directions (ms -1 ) 


nism (m) 

V 

specific volume of working gas (m 3 kg -1 ) 

L 

distance between two gears (m) 

V 

volume (m 3 ) 

Ldt 

length from the linkage 1 4 to top surface of the displac¬ 
er (m) 

w 

engine power (W) 

Lpt 

height from the linkage h to top surface of the piston 

Greeks 



(m) 

P 

density (kg m -3 ) 

m 

mass of working gas (kg) 

6 

wall thickness of the core cylinder (m) 

p 

pressure (Pa or N m 2 ) 

n 

engine efficiency 

q 

heat flux or heat transfer rate per unit area (W m -2 ) 

0 

crank angle 

Qin 

Qout 

input heat transfer rate (W) 
output heat transfer rate (W) 

CD 

engine speed (rad s -1 ) 

r, x 

two-dimensional cylindrical coordinates 

Subscripts 

f 

s 


u 

radius of the displacer (m) 


r 2 

R 

Re 

core radius of cylinder (m) 
gas constant (J kg -1 K -1 ) 

offset distance from the crank to the center of gear (m) 

solid 


features such as the geometries and dimensions of displacer and 
displacer cylinder, regeneration chamber, or the internal gas 
circuits, all of which can only be defined in multi dimensions. This 
means that they are not very useful in terms of the optimization of 
the geometries and dimensions of individual components during 
the design of a Stirling engine. In addition, heat transfer rates 
inside the expansion and compression chambers are often esti¬ 
mated by using constant convective heat transfer coefficients, a 
practice that has over-simplified the complexity of the real heat 
transfer mechanism taking place inside these chambers. These fac¬ 
tors together with other model assumptions result in poor predic¬ 
tive accuracy returned by these models. They tend to over-predict 
engine’s output power and efficiency by the margins from 30% to 
100%, or even more [8], 

A Stirling cycle involves very complicated heat and mass trans¬ 
fer processes, and the engine itself consists of many multi-dimen¬ 
sional components whose geometrical effects can be important. 
None of these can be resolved in great detail or very accurately 
by those aforementioned numerical models. Yet, resolving these 
complicated physical processes and the effects from some geomet¬ 
rical features is crucial to the prediction of the overall engine per¬ 
formance. Especially understanding the effects of the geometrical 
features is of practical importance because to design a new Stirling 
engine, engineers have to determine the detail geometries and 
dimensions of all components in the engine. Therefore, more 
sophisticated models are needed to improve the accuracy of the 
numerical analysis on Stirling engines. Computational fluid 
dynamics (CFD) is one of such numerical approaches by far that 
are capable of simulating those complicated processes taking place 
in Stirling engine and hence returning accurate prediction on the 
overall engine performance. It is also capable of resolving the 
effects caused by changes in small geometrical features and poten¬ 
tially be a great tool for the design and development of a new 
Stirling engine. However, there exist many challenges to apply a 
CFD analysis on a Stirling engine cycle. A Stirling engine cycle is 
unsteady and multi-scale in nature. In addition, the engine 
involves moving parts such as displacer and power piston; hence 
the volume of computational domain is constantly changing with 


time. These bring challenges in coding and physical modeling. In 
the aspect of coding, simulating the reciprocating motions of dis¬ 
placer and power piston requires the use of moving mesh tech¬ 
niques which involve complex algorithms. In the aspect of 
physical modeling, the gas flow is compressible; the heat transfer 
mechanism is conjugate heat transfer; and in most cases, the flow 
is turbulent. The combination of these factors gives rise to very 
complicated heat and mass transfer processes. To resolve these 
processes, it is necessary to use very fine grids and small time steps 
and solve many coupled governing equations at the same time, 
resulting in difficulty of convergence and long CPU time. Due to 
these challenges, many attempts to apply CFD analysis on Stirling 
engines focused on individual components such as the combustion 
chamber of the heater [12], So far, there have been few reports on 
using CFD approach to study the full cycle of a Stirling engine. 
Mahkamov [8] used both second-order approach and axisymmet- 
ric CFD approach to analyze the working process of a solar Stirling 
engine. They found significant differences in temperature data and 
engine performance obtained by both methods. The transient tem¬ 
perature variations by CFD approach were far from the harmonic 
distributions given by the second-order approach, and CFD model 
returned much accurate output power than the second-order 
model. Mahkamov [13] also reported applying a 3D CFD approach 
on a prototype biomass Stirling engine to improve its performance. 
The engine was also analyzed by the first-order and second-order 
models. The CFD analysis identified some key factors that reduce 
the power of the engine: a geometrical feature that causes high 
level of hydraulic losses in the regenerator, and the entrapment 
of the gas in the pipe connecting two parts of the compression 
space and the dead volume introduced by the pipe. Such detailed 
analysis on the effects imposed by multi-dimensional geometrical 
features can only be resolved by using a CFD approach. The 
research also highlighted the usefulness of adopting a CFD 
approach to identify and solve problems during the design of a Stir¬ 
ling engine. Zhigang et al. [14] used CFD approach to investigate 
the effect of a regenerator filled with porous-sheets heat exchang¬ 
ers. The CFD results were found in very good agreement with 
experimental data, further proving the capability of the CFD 
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approach to return accurate prediction on Stirling engine perfor¬ 
mance. Chen et al. [15] conducted a CFD study on the heat transfer 
characteristics of the working cycle of a 3D y-type Stirling engine. 
It was found that impingement is the major heat transfer 
mechanism in the expansion and compression chambers, and the 
temperature distribution is three dimensional and highly non- 
uniform across the engine volume at any given moment. 

Although CFD approach has been applied to the research and 
development of industrial Stirling engines, in open literature, there 
has not been a detailed CFD study on the basic heat transfer pro¬ 
cesses in the Stirling engine cycle. Yet, understanding the funda¬ 
mental processes of the Stirling engine cycle in great detail is a 
key step to improve the engine’s performance. The objective of this 
study is to perform CFD analysis on the working cycle of a /1-type 
Stirling engine to reveal the characteristics of its heat transfer 
processes. This particular type of Stirling engine cycle was chosen 
because of its large cycle power, high thermal efficiency, and 
popularity. It has been widely studied before by zero- or 
one-dimensional models, therefore comparison can be made 
between the results obtained in this study and those in literature. 


2. Mathematical model 


The Stirling engine analyzed in this study is the rhombic-drive 
/1-type Stirling engine developed by Cheng and Yu [16], where 
the engine was analyzed by an improved variant of the second- 
order model. The schematic of the geometrical configuration of 
the engine is given in Fig. 1. According to Fig. 1, the displacements 
of the piston and displacer are: 


Y.,( l ) — L pt + R d sin 6 -t- ^ 


L -- l l-R d cose 


( 1 ) 


x„(t) = L dt + R d Sine - 


It- [^- l j-R<,cos8 


and the velocities of piston and displacer are: 


u p (t) = ftdCiK cose - sine 


l\ - f 4 ^ - R d cos 0 


2 2 


xg-|-R d c °se 


')}• 


( 2 ) 


(3) 


u d (t) = R d co< cose + sine 


ll-(^- l j-R d cosd 


x(^- l j-R d cos e 


')}■ 


(4) 


The configuration of the computational domain is illustrated in 
Fig. 2. As seen, the domain is axisymmetric, allowing the engine 
cycle to be simulated by using two-dimensional mathematical 
model. The engine is very small with a core cylinder radius of just 
0.0205 m, thus the flow inside the engine can remain laminar even 
at relatively high engine speed. In addition, due to its small size, 
the engine is not equipped with a regenerator, and regeneration 
is performed by the working gas flowing forwards and backwards 
between the expansion and compression chambers through the 
narrow gap formed by the displacer and the inner wall of core cyl¬ 
inder: that is, the solid material of the core cylinder is acting like 
regenerator material to store and release heat. Such a simple 


Expansion 

chamber 



Fig. 1. The configuration and definition of geometric parameters of the /J-type 
Stirling engine. 
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Fig. 2. Schematic of the computational domain. 
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engine is ideal for the study of basic Stirling engine cycle. The 
assumptions made to facilitate CFD simulation are as follows: fluid 
viscosity and thermal properties of all materials are assumed 
constant, the working gas is air and assumed to be ideal gas, 
and effects from mechanical friction and thermal radiation are 
ignored. Therefore, the flow and heat transfer phenomena of 
the engine cycle can be solved by transient two-dimensional 
axisymmetric laminar and compressible Navier-Stokes equa¬ 
tions plus energy equation and ideal gas equation of state as 
follows [17]: 

Continuity equation: 

dp, d . 1 d 

l + ^ + r¥r^ = °- (5) 


Momentum equations: 


d 


1 d 


dt y 


( ^ U) + ax ( ^ U) + r Wr {p f r ' VU) 


r dr ' 


dp d fdu 
= Pg -]fx + ^{d-x 


u d { du\ p d ... 

+ 777 r- +L ( V .\7), 


r dr \ dr) 3 dx 


( 6 ) 


^Pf^+^PfUV)- 


19, - , 

-rdr^ rVV) 


dp ^(dv\ 
dr r dx\dx) 


u d 
r dr 


dv\ 

r dr) 


Energy equation: 

In gas: 

d d 1 d 

-(P / r / ) + -(p / t 1 r /)+? -(p / r,T» 


l 


dp 


c p/ L^t +v ' (pV)-p(v ' V) . 


+ '± 

c pf 


(V) 


d_fdTi 

dx\dx 


1 d 

+ --r- 


dT, 


r dry dr 



(dv 9u\ 2 

V9x + 9 rJ 



u = v= 0 - k fw =ks w Tf = Ts - (15) 

The outer wall of the core cylinder (at r=r 2 + 5) is assumed to 
maintain at a fixed temperature distribution: 

T s = T w (x). (16) 

3. Numerical Procedure 

The numerical procedure in this study is based on an in-house 
unstructured-mesh, fully collocated, finite-volume code, 
‘USTREAM’ developed by the corresponding author. This is the 
descendent of the structured-mesh, multi-block code of ‘STREAM’ 
[18], In a Stirling engine, the actions of the displacer and the piston 
cause expansion and compression of the engine volume, thus a 
moving mesh technique is needed to simulate such variations in 
volume. There are two commonly used moving mesh techniques: 
one is adding or deleting cells to change the volume of the compu¬ 
tational domain, and the other is to directly change the geometries 
of cells to expand or compress their volumes. In this study, the sec¬ 
ond strategy is adopted. The number of cells in the computational 
domain is fixed, and the expansion or compression of engine 
volume is performed by moving the local boundaries of cells to 
increase or decrease the volumes of the cells, and collectively, they 
can simulate the expansion or compression of the engine volume 
just like the action of an accordion. The motion of cell boundary 
creates local moving frame velocity components u and v. However, 
the component v is zero due to the fact that the motions of displac¬ 
er and power piston are both along the x-direction (axial direction). 

Another important issue is to simulate the compressible flow 
inside the engine. Since USTREAM is a pressure-based SIMPLE 
[19 finite volume code, the pressure correction equation has to 
be modified to account for variable density in compressible flows. 
The pressure-velocity coupling scheme in Lebon et al. [20] is 
employed here to perform this task. 


In solid: 


8 

at 


(PsTs) 


■a 

fdTA 

1 a 


dx 

ydxj 

+ 7dr 

V dr )J 


Equation of state: 


pv = RT f . 


( 8 ) 

( 9 ) 

( 10 ) 


where ii = u-u c , v = v - v c are relative velocity components 
between fluid and local moving frame that moves with u c and v c 
respectively in x and r directions. The last term on the right hand 
side of Eq. (8) accounts for viscous dissipation effect. As argued in 
Mahkamov [13], this effect is insignificant compared with other 
heat transfer effects in Stirling engine, thus it can be neglected. 


The initial conditions are: 

9 = 0°, u = v = 0, p = 101 kPa, T f — T s = T L . (11) 

Boundary conditions are:At x =X P : 
u = u p , v = 0, T f = T l . (12) 

At x = L t : 

u=v = 0,T f = T H . (13) 

The surface of displacer is assumed adiabatic, thus the condi¬ 
tions on displacer surfaces are: 

u = u d , v = 0. 5f = 0. (14) 


where n is the direction normal to the wall of displacer. On the 
solid/fluid interface (at r= r 2 ): 


4. Results and discussion 

The dimensions of the simulated /1-type Stirling engine are 
given in Table 1. The working gas is air with the gas constant of 
287 J kg~’ I< ’, and the core cylinder material is stainless steel with 
density, thermal conductivity, and specific heat of 7840 kg m 3 , 
43 W nr 1 K \ and 0.45 kj kg~’ K \ respectively. The initial crank 
angle and gas pressure are zero degree and 101 kPa, respectively. 
The external wall temperature distribution defined by Eq. (16) is 
assumed as: 

J 300 K, x < 0.08 m 1 

T(X) = | 300 K + ow^osX (800-300 K), x > 0.08 m J (1?) 


Table 1 

The dimensions of the geometrical parameters of the 
/?-type Stirling engine. 


r, (m) 

0.02000 

r 2 (m) 

0.02050 

G(m) 

0.00050 

1, (=l2-l3 = U)(m) 

0.01800 

k (m) 

0.15800 

L (m) 

0.04200 

i p[ (m) 

0.05093 

L dt (m) 

0.16374 

6 (m) 

0.07946 

R d (m) 

0.00360 

S (m) 

0.00100 
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Before the CFD simulation can be performed, the correctness of 
the code has to be validated. For this purpose, the engine was run 
under adiabatic condition at a very low engine speed of 60 rpm. 
The exact p-V relation in a reversible adiabatic cycle of air is 
pV 1 - 4 = C, where C is a constant. Since the current CFD code does 
not consider mechanical friction or viscous heat dissipation, the 
adiabatic cycle at low engine speed would be close to an ideal 
reversible adiabatic cycle. Fig. 3 shows the comparison of adiabatic 
p-V relations between the exact values and the present CFD 
results. It can be seen that the simulated p-V relation is almost 
identical to the exact relation, thus validating the correctness of 
the code. The next step of the CFD simulation is to determine a 
proper grid and time step interval to obtain grid and time-step 
independent solutions. In Cheng’s study, their baseline engine 
operated at an engine speed of 2000 rpm with hot and cold ends 
temperatures of 800 K and 300 K, respectively. An initial attempt 
to adopt the operation conditions of the baseline engine to perform 
the grid and time-step independency test was not successful. The 
solution converged up to the engine speed of 1860 rpm, and at 
2000 rpm, the code diverged at a certain time step. Further inves¬ 
tigation revealed that at 1860 rpm, the maximum velocity at the 
regenerative channel is about 44.0 ms~\ and the Reynolds num¬ 
ber, calculated by Re = >*£, is 1833, which is close to upper limit 
of laminar flow. Therefore, at the engine speed of 2000 rpm, the 
flow could have become turbulent within a certain period in an 
engine cycle, and the laminar flow model would fail to converge. 
For the purpose of grid and time-step independency test, a slightly 
slower engine speed of 1800 rpm was selected. Three different 
grids with cell numbers of 6264, 7864, and 9260, respectively 
were tested by using three different time-step intervals of 
8.333 x 10~ 5 s, 5.555 x 10~ 5 s, and 4.166 x 10~ 5 s, which respec¬ 
tively correspond to time-step numbers per cycle of 400, 500, 
and 800. Fig. 4 shows the grid with 6264 cells. Even though this 
is the coarsest grid, there are 20 cells allocated within the span 
of the regenerative channel. Cells are concentrated towards all 
solid boundaries to resolve the flow and thermal boundary layers. 
For each simulation, about 10 cycles were needed for the solution 
to become periodical, and each cycle took about 2-4 h in a PC, 
depending on grid size and time step number in a cycle. Judging 
from the values of maximum velocity, heat input, and work output 
within a cycle, grid and time-step independent solutions can be 
achieved by the grid with 7864 cells at 500 time steps per cycle, 



Fig. 3. Comparison of p-V relation between exact values and the present results in 
an adiabatic cycle with engine speed of 60 rpm. 



Fig. 4. Computational mesh. 



Fig. 5. Variations of total, expansion chamber, compression chamber volumes 
versus crank angle within an engine cycle. 


and this grid and the time step number have been used in this 
paper, in the following, the discussion will be focused on the case 
with the engine speed of 1800 rpm. 

Fig. 5 shows the variations of total, expansion chamber, and 
compression chamber volumes versus crank angle over an engine 
cycle. The maximum and minimum values of the total volume 
are respectively 3.090 x 10 5 m 3 and 1.774 x 10 5 m 3 , giving 
rise to a compression ratio of 1.742, and they occur at 0= 146.16° 
and 303.84°, respectively. It is also noticeable that the maximum 
volume of the compression chamber is larger than that of the 
expansion chamber, but the minimum volume of the former is 
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v (m 3 ) 


expansion chamber 



Fig. 6. Variations of averaged pressure versus volume in the expansion and 
compression chambers (a) current at 1800 rpm, and (b) Cheng and Yu [16] at 
2000 rpm. 


smaller than that of the latter. That is, the magnitude of volume 
variation in the compression chamber is larger. The p aV e~V dia¬ 
grams of the expansion and compression chambers are illustrated 
in Fig. 6(a). Here p ave is the averaged pressure in the individual 
chamber. Compared with the p-V diagrams in Cheng and Yu [16] 
(Fig. 6(b)), it has been found that apart from the differences in 
the magnitudes of minimum and maximum pressure, the tendency 
of p-V variation in the present study is very similar to that in 
Cheng’s study, which was calculated by an improved variant of 
the second-order model. This suggests that the simple second- 
order model is capable of predicting reasonably accurate pressure 
variations for the current engine. Since the pressure difference 
caused by gas flow inside the engine is just a tiny fraction of the 
total magnitude of pressure, the value of pressure is mainly 
determined by ideal gas equation (Eq. (10)). Therefore, the reason 
for this similarity should be due to the fact that pressure in the 
present CFD simulation and in Cheng’s second-order model is all 
obtained by the ideal gas equation. 

To demonstrate the general variations of temperature within 
the engine over an engine cycle, Fig. 7 shows the temperature con¬ 
tours of the entire engine at 10 different crank angles, from 0 = 36° 
to 360° with an increment of 36°. The motions of displacer and 
power piston pose different effects on the flow inside the engine. 
The combined motions of displacer and power piston cause gas 
to flow forwards and backwards between the expansion and 


compression chambers, while the motion of the power piston 
expands or compresses the total volume of the engine, which in 
turn introduces cooling or heating effects on the working gas. 
The gas flow that flows forwards and backwards between expan¬ 
sion and compression chambers introduces very strong heat con¬ 
vection effect in both chambers and results in complicated heat 
transfer behaviors between the gas and the solid boundaries and 
highly non-uniform temperature distributions across the chambers 
at any given moment. To this end, it is readily realized that the 
assumption of uniform temperature distribution and the adoption 
of constant convective heat transfer coefficients to account for the 
heat transfer effect in expansion and compression chambers in 
some second-order models is apparently too crude to represent 
the real heat transfer behaviors inside the engine. Therefore, the 
attention is directed towards the investigation on the engine’s heat 
transfer characteristics which can only be analyzed in detail by 
using CFD approach. 

From the expansion and compression chambers’ point of view, 
both go through one phase of receiving gas injection from the other 
chamber and the other phase of gas ejection to the other chamber 
through the regenerative channel. These phases are termed “injec¬ 
tion phase” and “ejection phase”, respectively. Basically, the injec¬ 
tion phase of the expansion chamber is the ejection phase of the 
compression chamber, and vice versa. In the following, the discus¬ 
sion on the flow and heat transfer characteristics is divided into 
injection phase, ejection phase, and regenerative channel 
subsections. 

4.J. Injection phase 

The injection phase in the individual chamber generally coin¬ 
cides with the expansion process of its chamber volume. In the 
expansion chamber, this is from 0 = 216° to 22° (or 382°); and in 
the compression chamber, it is from 0° to 198°. As shown in 
Fig. 7, this phase in the individual chamber can be characterized 
by the appearance of a major thermal plume in that chamber. 
However, there is a brief interval, from 0° to 22°, when gas at the 
lower part of regenerative channel is injected into the compression 
chamber, and gas at upper part is injected into the expansion 
chamber. This peculiar phenomenon is due to the downwards 
motion of displacer that continues to drag gas to the compression 
chamber and the onset of the expansion process of the expansion 
chamber that creates vacuum and sucks the gas into expansion 
chamber. The injection phase poses significant impact on the heat 
transfer in both chambers. To demonstrate this effect, Fig. 8 shows 
the enlarged views of velocity vectors and temperature contours as 
well as local heat flux distributions in the expansion chamber at 
6 = 298.8°. Fig. 8(a) shows that gas injection produces a plume of 
cold gas that impinges on the upper heated wall of the chamber 
(the hot-end wall). This impingement brings the cold gas in direct 
contact with the hot-end wall, resulting in impingement heat 
transfer which is well known to be the most effective convective 
heat transfer mechanism. Therefore, the cold gas get injected into 
the expansion chamber is rapidly heated up. This can be verified 
from the distributions of the local heat flux per unit area at 
6 = 298.8° along the surfaces of power piston, regenerative channel, 
and hot end, Si, s 2 , and s 3 , given in Fig. 8(b). Here, the surfaces Si, s 2 , 
and s 3 are defined in Fig. 2, and local heat flux q is calculated by: 


Under this definition, the heat transfer from solid wall (sur¬ 
rounding) to gas (thermodynamic system) is positive, and that 
from gas to wall is negative. This particular phase angle was chosen 
because it corresponds to the maximum local heat flux at the hot- 
end wall (s 3 ). Comparing Fig. 8(a) and (b), it can be seen that the 
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Fig. 7. Temperature contours across the whole engine at 10 different crank angles. 


maximum local heat flux (in Fig. 8(b)) occurs within the vicinity of 
the impingement point (in Fig. 8(a)) of the cold gas jet on the hot- 
end wall and its value reaches a phenomenal value of 80 kW m 2 , 
which is also the absolute maximum value during the entire 
engine cycle. As seen in Fig. 8(a), the jet of cold gas forms a primary 
vortex which in turn induces multiple secondary vortices at its 
sides. These secondary vortices promote mixing of cold and hot 
gas and further facilitate convective heat transfer in the expansion 
chamber. Their effects can be identified by the smaller peaks at 
both sides of the primary peak on s 3 ; however, their magnitudes 
are much smaller than the primary peak produced by the impinge¬ 
ment effect. Near the center of expansion chamber, local heat flux 
is very small as convection becomes very weak. The local heat flux 
distributions prove that the dominant heat transfer mechanism in 
the injection phase of the expansion chamber is the impingement 
of the cold gas on the hot-end wall. 

Fig. 9 shows the enlarged views of velocity vectors, temperature 
contours, and local heat flux distributions of the compression 


chamber at 9= 172.8°. Similar impingement phenomenon in the 
expansion chamber can also be observed during the injection 
phase of the compression chamber, but the impingement results 
in rapid cooling of the hot gas that enters the compression cham¬ 
ber. The maximum (in terms of absolute value) local heat flux 
along the cold-end surface Sj during the injection phase is 
-24.6 kW m~ 2 . This is much smaller than that in the expansion 
chamber during its injection phase. In addition, unlike what hap¬ 
pened in the expansion chamber, it is neither the maximum local 
heat flux on S! during the entire engine cycle. The real maximum 
local heat flux on s, occurs in the ejection phase which will be dis¬ 
cussed in the next subsection. 

4.2. Ejection phase 

The ejection of gas in expansion or compression chambers is 
mainly caused by the reduction of volume in the individual cham¬ 
ber. By referring to Fig. 5, this phase is from 0 = 53° to 215° in the 
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6=298.8° 



Fig. 8. Temperature contours, velocity vectors, and local heat flux distributions in 
the expansion chamber at 0 = 298.8°; (a) temperature contours and velocity vectors 
and (b) local heat flux distributions. 



Fig. 9. Temperature contours, velocity vectors, and local heat flux distributions in 
the compression chamber at 8 = 172.8°; (a) temperature contours and velocity 
vectors and (b) local heat flux distributions. 


expansion chamber and from 0= 199° to 359° in the compression 
chamber. As mention earlier, the reduction of chamber volume is 
due to the motions of the displacer and power piston. When 
displacer moves towards the hot end, the volume of expansion 
chamber is reduced. On the other hand, when displacer and power 
piston move toward each other, the volume of compression cham¬ 
ber is reduced. In this manner, the displacer and power piston 
reduce the volume in the individual chamber by decreasing the 
distance between the upper and lower boundaries of that chamber, 
and this is equivalent to pushing the gas against the solid walls. In 
the expansion chamber, the upper wall is the hot end, while in the 
compression chamber the lower wall is the cold end. The pushing 
of gas by the displacer and power piston creates an effect similar to 
impingement that largely compresses the isothermal lines towards 
the hot or cold ends and increases the local temperature gradient. 
As a result, it also promotes heat transfer in a way similar to flow 
impingement does. It can be seen in Fig. 5 that in the compression 
chamber, the ejection phase mostly coincides with the compres¬ 
sion process of the entire engine volume, and the negative heat 
transfer is further boosted by the compression effect which heats 
up the gas in the compression chamber and creates even larger 
temperature difference between the gas and the cold-end wall. 


The combined heat transfer effects can be demonstrated by veloc¬ 
ity vectors, temperature contours, and local heat flux distributions 
in the compression chamber at 0 = 230.4° shown in Fig. 10. The 
temperature contours in Fig. 10(a) show that the core temperature 
in the compression chamber reaches 430 K due the compression 
process, and the isothermal lines are very dense towards the 
cold-end wall due to the motion of power piston. Consequently, 
the local heat flux on S] is highly elevated as can be seen in 
Fig. 10(b). The maximum absolute value reaches 50 kW m~ 2 , which 
is more than twice the maximum value in the injection phase and 
is the absolute maximum value on Si. In the expansion chamber, on 
the other hand, the ejection phase does not coincide so well with 
the expansion process of the entire engine volume; therefore, the 
promotion on positive heat transfer by this mechanism is not very 
significant. On the contrary, Fig. 5 shows that from 6 = 150° to 220°, 
the ejection phase of the expansion chamber coincides with the 
compression of total engine volume, thus the gas inside the expan¬ 
sion chamber is heated not only by the hot-end wall but also by the 
compression effect. As a result, the maximum temperature in the 
expansion chamber reaches as high as 860 K, which is higher than 
the hot-end temperature. The elevated gas temperature then gives 
rise to negative heat flux on s 3 as can be seen in Fig. 9(b). This indi¬ 
cates that heat output is not confined to the cold end but can occur 
at hot end within a short interval during an engine cycle. 
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Fig. 10. Temperature contours, velocity vectors, and local heat flux distributions in 
the compression chamber at 0 = 230.4°; (a) temperature contours and velocity 
vectors and (b) local heat flux distributions. 


4.3. Regenerative channel 

The major function of a regenerative channel is to pre-heat the 
gas during the injection phase of the expansion chamber and to 
pre-cool the gas during the injection phase of the compression 
chamber. From the local heat flux distributions shown in 
Figs. 8(b) and 9(b), it can be seen that the regenerative channel per¬ 
forms this task relatively well as the heat flux on s 2 in Fig. 8(b) is 
mostly positive (heating gas), and that in Fig. 9(b) is mostly nega¬ 
tive (cooling gas). However, there are exceptions. In the ejection 
phase of the compression chamber, Fig. 8(b) shows that at the 
lower section of s 2 (near entrance of the regenerative channel), 
heat flux is negative. This happens when the gas with a tempera¬ 
ture higher than the regenerative-channel wall temperature from 
the compression chamber enters the regenerative channel. The 
higher gas temperature in the compression chamber is due to the 
heating of gas by the compression effect mentioned in Section 4.2. 
Exception also occurs in the ejection phase of the expansion cham¬ 
ber. Fig. 11 shows that distributions of local heat flux and velocity 
vectors and temperature contours in the vicinity of the entrance of 
the regenerative channel in the expansion chamber at 6 = 93.6°. In 
Fig. 11(a), a large spike of positive heat flux is visible downstream 
the entrance region, indicating that gas is heated in this part of 
regenerative channel. The reason is due to the discharge of some 
low temperature gas into the regenerative channel as can be seen 
in Fig. 11(b). The low-temperature gas exists because of the non- 
uniform heating of gas in the expansion chamber during its injec¬ 
tion phase. Although very effective, the impingement heat transfer 
mentioned in Section 4.1 is also highly concentrated at some 
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Fig. 11. Temperature contours, velocity vectors, and local heat flux distributions 
near the entrance of regenerative channel in the expansion chamber at 0 = 93.6°; (a) 
local heat flux distributions and (b) temperature contours and velocity vectors. 


localized region. As seen in Fig. 7, this results in a pocket of gas 
in the expansion chamber that has not been well heated and 
remains at relatively lower temperature. Fortunately, these 
exceptions only appear in small localized regions of the regenera¬ 
tive channel and they only last for some brief moments. The regen¬ 
erative channel actually functions as it is designed to function at 
the majority of time in an engine cycle. In the following, attention 
will be focused on discussing average and integrated quantities 
such as average temperature variation, heat input, and work 
output within an engine cycle. 

4.4. Average and integrated quantities 

Since some second-order model adopt a very simple practice to 
account for heat transfer, it is interesting to see what level of accu¬ 
racy such practice can achieve. Cheng and Yu [16] has illustrated 
the distributions of temperature and rates of heat transfer in the 
expansion and compression chambers. As mentioned earlier, those 
results are for engine speed of 2000 rpm, which is about 11% higher 
than the current engine speed of 1800 rpm. However, we can still 
compare the similarity in the tendency of the variations of some 
quantities returned by their second-order model and the current 
CFD approach. The distributions of average temperature and rate 
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of heat transfer in the expansion and compression chambers by the 
current study are given in Figs. 12 and 1 3, respectively. Results from 
Cheng’s study were also digitalized and plotted in these figures for 
comparison. Note that Cheng used time as the x-axis in their plots. 
Fig. 12 shows that there appears some similarity in the tendency of 
temperature variations. Yet, there are some major differences. One 
is the minimum temperature in the expansion chamber happens at 
different crank angle; and another is that the current study predicts 
a prolong period of low temperature that results in almost zero heat 
transfer in the compression chamber, which is not captured in 
Cheng’s results. Fig. 13 shows that the similarity in the tendency 
of heat-transfer-rate variation is much worse. Fig. 13(a) shows 
the presence of multi peaks and valleys in the heat-transfer-rate 
variations in both chambers within an engine cycle. The results 
from the second order model in Fig. 13(b), on the other hand, only 
predicts a single peak and a single valley in the expansion chamber 
and two peaks and two valleys in the compression chamber. The 
complex variations in heat transfer rate are rooted in the compli¬ 
cated heat transfer mechanisms that result in highly non-uniform 
temperature distributions mentioned in the earlier sections. Using 




Fig. 12. Variations of average temperature in the expansion and compression 
chambers; (a) current at 1800 rpm, and (b) Cheng and Yu [16] at 2000 rpm. 



Expansion 



Fig. 13. Variations of heat transfer rates in the expansion and compression 
chambers; (a) current at 1800 rpm, and (b) Cheng and Yu [16] at 2000 rpm. 


the data in Figs. 12 and 13, heat transfer coefficient can be evaluated 
by: 

h = T Q t , ( 19 ) 

* w * ave 

where T w is either the hot-end temperature T H or the low-end tem¬ 
perature T L respectively at the expansion or compression chamber. 
The distributions of heat transfer coefficient obtained are illustrated 
in Fig. 14. As seen, the heat transfer coefficients are not uniform in 
both chambers, and there are some anomalies of huge oscillations 
from positive to negative (or negative to positive). The anomalies 
indicate that simple equation such as (19) would fail to account 
for the heat transfer effect when applied to the CFD results. To this 
end, it is clear that the practice of using constant heat transfer coef¬ 
ficient to calculate heat transfer rate in some second-order models 
fails to account for the effects produced by impingement, mixing of 
gas, and non-uniform heating, and hence, such a model is not capa¬ 
ble of predicting rates of heat transfer very accurately in both 
chambers. 
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Fig. 14. Distributions of convective heat transfer coefficients calculated by Eq. (19). 

From the discussion in the previous subsections, it can be real¬ 
ized that heat input and output can happen at the hot end, cold 
end, and in the regenerative channel. Therefore, the integration 
of heat input cannot be confined just along the wall of the 
expansion chamber. Neither can the integration of heat output be 
confined along the wall of compression chamber. The rates of heat 
input and output are calculated by integration along all solid 
boundaries: 

r.s pt=t 0 +t p /*S=Si +S2+S3 

Q.in = T~ / 2 qdAdt, if q > 0. (20) 

Jt=t 0 JS=0 

/-.-i rt=to+tp rS=s-[ +S2+S3 

Q ouf = — / / 2 qdAdt, if q < 0, (21) 

where the upper and lower limits in the inner integration mean that 
the integration is performed along walls of compression chamber, 
regenerative channel, and expansion chamber. In terms of the eval¬ 
uation of output power, the focus is on the surface of the power pis¬ 
ton because this is the exact place engine work is delivered. The 
power is calculated by: 

W = — <j> j p2nrdrdx, (22) 

where the negative sign on the right hand side of equation is to 
account for the fact that positive work occurs when dx is negative, 
and the circular integration is from minimum X p to maximum X p 
and back to minimum X p . The calculated rates of heat input, heat 
output, and power are 82.5 W, 69.7 W, and 11.7 W, respectively, 
resulting in a thermal efficiency q of 14.2%. Here, q is calculated 
simply by dividing the engine power by the rate of heat input. Ide¬ 
ally, the total heat flow in the regenerative channel should be zero. 
The heat flow rate along the wall of the regenerative channel in this 
case is calculated as 1.512 W, which is very small compared with 
the rates of heat input. 

To help readers understand the dynamic processes inside the 
current Stirling engine, five videos have been made. The files 
named “temperature”, “vector”, “flux”, “expansion”, and 
“compression” respective play the transient variations of tempera¬ 
ture contours, velocity vectors, local heat flux, and enlarged 
temperature contours and velocity vectors in expansion and com¬ 
pression chambers within an engine cycle. 


5. Conclusion remarks 

An in-house CFD code has been developed and applied to a 
rhombic-drive p- type Stirling engine to study the heat transfer 
characteristics in the thermal cycle of the Stirling engine. The 
CFD approach allows the heat transfer characteristics to be inves¬ 
tigated in great detail. The results include temperature contours, 
velocity vectors, and distributions of local heat flux at several 
important time steps in an engine cycle as well as average temper¬ 
ature variations, integrated rates of heat input, heat output, and 
engine power. Some conclusion remarks from this study can be 
drawn as follows: 

1. The complicated heat transfer behaviors in the engine result in 
highly non-uniform temperature distributions in all chambers 
at most of the time in an engine cycle. The assumption of uni¬ 
form chamber temperatures and the use of constant heat trans¬ 
fer coefficients to account for heat transfer adopted by some 
zero- or one-dimensional models are too crude to reflect the 
reality. 

2. In the injection phase, impingement heat transfer is the major 
heat transfer mechanism which allows gas to be heated (in 
expansion chamber) or cooled (in the compression chamber) 
very effectively. However, this effect tends to concentrated at 
some localized region, resulting in non-uniform heating or cool¬ 
ing of gas. 

3. In the ejection phase, the motions of displacer and power piston 
also create a sort of heat transfer effect similar to impingement, 
but this effect acts on a much larger region and also promotes 
heat transfer effectively especially in the compression chamber 
where the negative heat transfer is further boosted by the com¬ 
pression of the entire engine volume. 

4. The regenerative channel functions as it was designed for at 
most of the time within an engine cycle. It pre-heats or pre¬ 
cools gas as gas is injected into the expansion or compression 
chamber. However, exceptions can occur near both regenera¬ 
tive channel entrance regions at some brief moments. 

5. Rates of heat transfer in both expansion and compression 
chambers cannot be predicted very accurately by those sec¬ 
ond-order models which adopt a very simple practice to 
account for heat transfer effect. 

The future work includes: 1. a parametric study on some 
important geometrical factors to understand their impact on 
engine performance and 2. the introduction of a regenerator filled 
with porous medium into the engine model to more realistically 
account for the regeneration effect. 
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